The Spin-1 Heisenberg Antiferromagnet: New Results from Series Expansions 
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We calculate ground state properties (energy, magnetization, susceptibility) and one-particle spec- 
tra for the S = 1 Heisenberg antiferromagnet with easy-axis or easy-plane single site anisotropy, 
on the square lattice. Series expansions are used, in each of three phases of the system, to obtain 
systematic and accurate results. The location of the quantum phase transition in the easy-plane 
sector is determined. The results are compared with spin-wave theory. 
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I. INTRODUCTION 

Magnetic materials with S = 1 ions have been of inter- 
est for many years. The 'classic' 2-dimensional Heisen- 
berg antiferromagnet K2NiF4 was studied in the 1970s 
[l|. In the 1980s a number of (weakly coupled) linear 
chain systems were investigated, including CsNiCla [2|, 
which has a weak axial anisotropy, CsFeBrs [J], which 
has strong planar anisotropy, and the complex mate- 
rials NENP (Ni(C2H8N2)2N02(C104)) Q and NENC 
(Ni(C2H8N2)2Ni(CN4)) ;5|, which have weak and strong 
planar anisotropy, respectively. The spin gaps observed 
in the weakly anisotropic materials [2, 4] are believed to 
be examples of the behaviour predicted by Haldane ia|. 
More recent work includes molecular oxygen adsorbed on 
graphite [3], the bilayer material Ba3Mn2 08 @,@], a new 
spin gapp ed material NiGa2S4 [Ic|, which, it has been ar- 
gued [IllliJl, may be a 'spin nematic' Ij|, and a system 



of spin-1 bosonic atoms in an optical lattice [1^ 

We have previously used series expansion methods to 
study a wide range of spin-1/2 quantum antiferromag- 
nets [15'|. Quantum fluctuations will be reduced when 
S = 1, but new physical features are possible. Besides 
the gapped Haldane phase in 1-dimensional systems p, 
additional terms, such as biquadratic exchange and/or 
single site anisotropy, which are absent in spin-1/2 sys- 
tems, can lead to quadrupolar/nematic phases with long 
range order but no magnetic moment, novel quantum 
phase transitions [1^, and richer excitation spectra. We 
explore some of these issues in the present work, within 
the context of the Hamiltonian 



H = jJ2s,-Si-Dj2{S! 



(1) 



<ij> 



which describes a Heisenberg antiferromagnet (J > 0) 
with isotropic exchange and a single ion anisotropy which 
gives rise to an easy axis {D > 0) or easy plane (£> < 0). 
A significant single ion anisotropy is believed to be 
present in many of the 1-dimensional materials referred 
to above, and has been included in analyses of the exper- 
imental data. Consequently there has been much theo- 
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FIG. 1: The phase diagram of the spin-1 J-D model on the 
square lattice at zero temperature, showing the Ising anti- 
ferromagnetic phase (lAFM), planar antiferromagnetic phase 
(PAFM), and quantum paramagnetic phase (QPM). 



retical work devoted to the Hamiltonian ([1]) on a linear 
chain [l7|, llM UM- ^^^' higher dimensions various ap- 
proaches have been used, including mean-field type the- 
ories 2(i, i21j, spin wave approximations |22l . |23[ and a 
coupled cluster calculation [2j|. We will compare our 
results with each of these, where possible. 

The present work deals with the 2-diniensional square 
(SQ) lattice. For D = the spin-1 system has Neel or- 
der at zero temperature [25|, with quantum fluctuations 
reducing the staggered magnetization by some 20% [2^ . 
When D > the system will order along z, the 'easy 
axis'. The order parameter will be 'Ising-like' and long- 
range order will persist at finite temperature, up to a crit- 
ical line Tc{D) with Ising (n=l) exponents. For D < 0, 
on the other hand, the z-axis is a 'hard' direction and 
the spins will order antiferromagnetically along some di- 
rection in the x-y plane (at least for small |D|). The 
order parameter has n=2 components and the continu- 
ous symmetry will be spontaneously broken. Long-range 
magnetic order will not persist to finite temperature (the 
Mermin- Wagner theorem) although one may expect a 
Kosterlitz-Thouless transition. For large negative D the 
physics will be quite different. In the limit D — + — cx3 the 
ground state will be a simple product state with S^ — 
at all sites. This is a quadrupole state with no mag- 
netic order. Thus we expect a quantum phase transition 
at some D = Dc, which we will locate using our series 
approach. The phase diagram is illustrated in Figure [TJ 

It is also of interest to study the elementary excitations 
above the ground state. For D = these are magnons, 
with Goldstone modes at k = (0,0) and (7r,7r). A gap 



will open in the spectrum for small positive D, which 
will be proportional to \fD, according to spin-wave the- 
ory. On the other hand for large \D\ the picture will 
be quite different. For large negative D the excitations 
will consist of isolated S*^ = ±1 spins, which have been 
termed excitons and anti-excitons [18|. These will have 
an energy gap, which we expect will vanish as Z? — *■ Dc— • 
For large positive D an excitation with A5^ = 2 (i.e. 
5^ = — 1^^5^ = -|-1) will have a smaller energy than 
a single magnon. We expect, and confirm below, that 
there is a 2-magnon bound state, and we calculate its 
dispersion curve. 

To conclude this Introduction we will briefly describe 
the series expansion approach, for both ground state bulk 
properties and excitations, referring the reader to our 
recent book Il5| for further details. The approach is 
based on writing the Hamiltonian in the usual pertur- 
bative form 



ii = H. 



\V 



(2) 



where Hq has a simple known ground state. This subdi- 
vision is carried out in various ways, appropriate to the 
different phases of the model. Series are derived for the 
ground state energy, order parameter and other quanti- 
ties of interest, in powers of A, and extrapolated to A = 1 
by numerical methods (Fade approximants or integrated 
differential approximants). Calculations are carried out 
for a sequence of finite connected clusters, and the results 
combined to obtain bulk lattice properties. 

A similar approach is used for excitations. An orthog- 
onal transformation is used to obtain an 'effective Hamil- 
tonian' matrix for each cluster, yielding transition ampli- 
tudes for the excitation in real space, which are then com- 
bined to obtain dispersion curves throughout the Bril- 
louin zone. Foints where the gap approaches zero can be 
easily identified, and corresponding series obtained for 
the gap itself. 

In the body of the paper we will present and analyse 
various results in the easy-axis (I) > 0) phase (Section 
II), and in the easy-plane {D < 0) phase (Section III). 
Section IV contains a summary and conclusions. 



II. THE EASY AXIS {D > 0) CASE 
A. Bulk Ground State Properties 

We have computed series for the ground state energy 
per spin and for the ground state staggered magnetiza- 
tion to order A^^, where Hq and V are 

Ho = JY. ^^^j - DY^^Stf ~ h J2 mSl (3) 

< ij > i i 

i^(s+5- + sr^;) (4) 
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FIG. 2: Ground state energy per spin for the easy-axis J- 
D model on the SQ lattice. The individual points are the 
estimates from series . The various lines are results of different 
spin-wave approximations: SWl, short-dashed line; SWla, 
long-dashed line; SW2, solid line. For convenience, we plot 
Eo + D versus D, and set J — 1. 



for various values of D. Here /i is a staggered field, in- 
cluded to allow calculation of the order parameter, and 
rji = ±1 on the respective sublattices. We do not present 
the series coefBcients here, but can provide them on re- 
quest. 

Figures [Hand [3] show the ground state energy and stag- 
gered magnetization versus D. 

For purposes of comparison, we also present results 
from conventional spin- wave theory, first order spin wave 
theory (SWl), modified first-order theory (SWla) and 
second order theory (SW2). The first order theory gives 



Eo = -4J-(l + 77)i^ + l^ 
3 1 ^ 

9 Iv 2^ 



M, = - - 

2 N 



N ■ 
4 J + 2r]D 

Cfe 



£k 



where e^, the spin- wave energy, is given by 



Ck = U^{l + iiD/{2J)f-^l 



(5) 
(6) 

(7) 



V 



<tj> 



with 7fc — (cosfcj, -|- cosky)/2. The constant 77 — 1,1/2 
for SWl, SWla respectively. This factor arises from a 
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FIG. 3: Staggered magnetization for the easy-axis J-D model 
on the SQ lattice. Notation as in Figure [S] 



FIG. 4: Single-magnon excitation energy at D/J = 1.0 along 
symmetry lines in the Brillouin zone, obtained from series (full 
dots) and spin wave approximations. Notation as in Figure 



choice in treating the anisotropy term. In large S theory 
77 = 1 but for S = 1, as here, normal ordering of the 
quartic boson terms gives rj = 1/2. This is explained in 
Appendix A, where we also present the more complex 
second-order theory. 

As is apparent from Figure [21 SWl is a rather poor 
approximation, but SWla and SW2 are in almost quan- 
titative agreement with the series results, within 1%. For 
the magnetization, a simlar conclusion applies. Note that 
the magnetization exhibits a square-root cusp behaviour 
near D = 0, in agreement with the spin- wave theory. 



B. 1-Magnon Excitations 

At least for small D the lowest energy excitations, in 
the unperturbed system, consist of a single spin excited 
from its ordered S^ = ±1 state to S^ = 0, i.e. AS^ = ±1. 
The quantum fluctuations, embodied in the perturbation 
V, will then allow this to propagate through the lattice, 
forming a coherent magnon band. 

Within spin wave theory the excitation energy is given 
by equation ([7]) (or the more general result in Appendix 
A). We have computed a series expansion for the excita- 
tion energy to order A^^, following the original work of 
Gelfand [27| (see also ref. 13). 

In Figure [3] we show the excitation energy, along sym- 



metry lines in the Brillouin zone, for the case D/J = 1, 
and again for comparison the spin- wave results. Again 
we see that SWl is a rather poor approximation, but 
SWla and SW2 provide an excellent description of the 
data, with SWla actually better than SW2. 

The dispersion curves are smooth and rather feature- 
less. The most significant feature is the opening of a gap 
at k = (0, 0) for any non-zero D. This reflects, in the 
easy axis case, the fact that the remnant 0(2) symmetry 
of the Hamiltonian is not spontaneously broken in this 
case, and so Goldstone modes are absent. 

Figure [5] shows the dependence of the gap at k = (0, 0) 
on the coupling D. The VD dependence at small D pre- 
dicted by spin wave theory is clearly evident. At large 
D the single-magnon gap is predicted to increase linearly 
with D. 



C. Large D Excitations 

It is clear that for large D the single-magnon excita- 
tions of the previous subsection will not be the lowest 
energy excitations of the system. Their energy will be of 
order D, whereas an excitation with AS*^ = 2 will have 
an energy of order J. Such an excitation, created at a 
particular site, can again propagate through the lattice. 
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FIG. 5: Single-magnon excitation energy at k = (0, 0) as 
a function of D in the easy-axis region (setting J — 1). The 
solid line is the second order spin wave approximation (SW2). 



FIG. 6: Dispersion relation of the AS^ = 2 excitation at 
D/J = 1.0. The points with error bars are the series esti- 
mates. The dashed line is the lower limit of the 2-magnon 
continuum. 



forming a quasiparticle band. We may think of this as a 
two-magnon bound state where the magnons are bound 
on the same site (see Appendix A). 

Figure [HI shows the dispersion relation for the A5^ = 2 
excitation at D/J = 1 compared with the lower edge 
of the two-magnon continuum. It can be seen that in 
the mid-region of the plot the excitation indeed seems 
to lie below the continuum, becoming a bound state at 
slightly below this coupling. The energy here is close 
to the asymptotic limit of 8 units. In the wings of the 
plot the error bars are much larger, and the excitation 
may not be bound: it is possible that these facts may be 
related, At higher values of D/J, the bound state energy 
remains close to 8 units, so that the binding energy rises 
almost linearly with D/J. 



Xs, for various values of D. These are then analysed via 
standard Dlog Fade approximants to obtain the critical 
temperature and exponent. Figure [7] shows the critical 
temperature versus D and, for comparison, the mean field 
result [21j. One would not expect MFA to give accurate 
results in 2 dimensions, and indeed there is a sizeable 
discrepancy. We note that MFA gives a finite critical 
temperature even for the isotropic case, D = 0, which 
violates the Mermin- Wagner theorem. The critical expo- 
nent 7 is consistent with the universal 2D Ising exponent 
7/4, although there is substantial scatter in the estimates 
from these relatively short series. 



III. THE EASY-PLANE (D < 0) CASE 



D. Finite Temperature Phase Transition 

Since the continuous 0(3) symmetry of the Heisen- 
berg model is destroyed by an easy-axis anisotropy term, 
and the order parameter is Ising-like, taking one of two 
possible values, the ordered ground state will persist to 
finite temperature, up to a critical temperature Tc{D). 
We have derived high-temperature series in the variable 
K = J /ksT to order K^"^ for the staggered susceptibility 



The easy-plane case shows much more interesting 
physics, including, as we shall see, two distinct phases 
separated by a quantum phase transition (QPT). 

For small \D\ the spins will be preferentially in the 
x-y plane (choosing z as the hard axis) and the Hamilto- 
nian will have 0(2) symmetry. At T = this symmetry 
will be spontaneously broken and the system will exhibit 
Neel order in some direction, reduced by quantum fluctu- 
ations. We refer to this as the planar antiferromagnetic 
phase (PAFM). The broken 0(2) symmetry will result in 
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FIG. 7: Critical temperature versus D j J for the S = 1 easy- 
axis model on the SQ lattice. The filled circles are the series 
results, with estimated errors no larger than the size of the 
symbols. The line is a guide to the eye. The dashed line is 
the MFA result [2H. 



FIG. 8: Ground state energy per site in the easy-plane region 
for the S = \ J-D model on the SQ lattice. The squares are 
series results for the PAFM phase, and the the dashed line is 
the first order spin wave estimate in the PAFM phase. We 
have set J = 1. 



a single gapless Goldstone mode. In the following we will 
present results from series expansions for both ground 
state bulk properties and for single-magnon excitations. 
These will again be compared with spin wave theory. Al- 
though there will be no ordered phase at finite tempera- 
ture we expect a finite temperature Kosterlitz-Thouless 
transition. However we do not explore this aspect. 

For large \D\i where the anisotropy term is dominant, 
we expect the system to prefer a singlet phase where spins 
are in the S*^ = state. This phase has no magnetic or- 
der, and is aptly referred to as a quantum paramagnetic 
phase (QPM). Quantum fluctuations, arising from the 
exchange terms, will modify the state. Low energy exci- 
tations in the QPM phase consist of spins excited to the 
5^ = ±1 states, which have been termed 'excitons' and 
'anti-excitons'. In the following we derive series for both 
bulk properties and excitations in the QPM phase. An 
analytic approximation, due to Papanicolaou [l8| , is also 
presented and compared with the series results. 

As \D\ is reduced in the QPM phase (or increased in 
the PAFM phase) a quantum phase transition is found 
to occur and we use series expansions to locate this tran- 
sition accurately, and to study its properties. 



A. PAFM Phase: Bulk Ground State Properties 

It is convenient to rotate the spin axes and to write 
the Hamiltonian as 



-:x\l 



= JY. ^^^0 + ^^ E ^^tSj + SrS+) 

-\DYisr+sn' (8) 



where z is the ordering axis and x is the hard axis. We 
now decompose H — Hq + XV, with 



Hn = 



V = 



J Y sts; + \d J2is-r h Y msi (9) 



<u> 



\DY{StSt + SrSr) 



(10) 
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FIG. 9: Staggered magnetization in the PAFM phase of the 
5 = 1 J-D model on the SQ lattice. The points are series 
estimates, and the short-dashed line is first order spin wave 
theory. The long-dashed line is a fit in the critical region. 



where we have dropped a constant term —ND and, as 
before, included a staggered field term. 

We have computed series for the ground state energy 
and staggered magnetization to order A^^. Figures [8] and 
inishow the ground state energy and magnetization versus 
D, as obtained by analysis of these series. 

For comparison we show the results of first order spin 
wave theory (SWl). The details of this are given in Ap- 
pendix B. In Figure [3 we see that the SWl theory gives 
quite a good description of the bulk ground state energy 
in the PAFM phase, as compared with the series esti- 
mates (squares). 

The series estimates for the magnetization in the 
PAFM phase were obtained as follows. The perturbation 
series in A typically exhibit a singularity of form (Ac — A)*^ 
for Ac only a little larger than 1 , so a naive Pade extrapo- 
lation in A gives unreliable results. Instead, at each value 
of D/J we estimate Ac and a using standard Dlog Pade 
methods, and then extrapolate the series to A = 1 using 
a variable 5 = 1 — (1 — A/Ac)"^. The data seem to show 
a crossover from a singularity with very small a at small 
negative D, to another one with a ~ 0.27 nearer the crit- 
ical point. The estimated magnetization at the crossover 
point, around D/J = —3.5, shows large error bars. 

The results are shown in Figure [51 We note that SWl 



FIG. 10: Single-magnon dispersion curves for the PAFM 
phase for D/J = —1.0 (circles) and D/J = —5.0 (squares) 
The dashed fines are from the first order spin wave theory. 



gives qualitatively the right behaviour at small \D\, but 
becomes rather poor for large \D\, where it shows no sign 
of the decrease towards the transition point. Near Z? = 0, 
the magnetization shows a square-root cusp, the mirror of 
that in the easy-axis region, marking the transition from 
the easy-axis to the easy-plane phase. The series results 
show a rapid decrease in magnetization beyond D/J ^ 
—5, heralding the expected quantum phase transition to 
the QPM phase. However the error bars are large and it 
is not possible to locate the transition with any degree of 
precision from the magnetization alone. The fit in this 
region will be discussed in Section IIIIDI 



B. PAFM Phase: Elementary Excitations 

We have computed series for the single-magnon exci- 
tations in the PAFM phase, to order A^°. The series 
are analysed to compute the magnon energies e(k), and 
these are shown in Figure [TOl along symmetry lines in the 
Brillouin zone, for values D/J — —1.0 and —5.0. Again, 
for comparison, we show the result of first order spin 
wave theory (Appendix B). We note that spin fluctua- 
tions transverse to the ordering direction are no longer 
isotropic, and hence the full Brillouin zone must be used. 

The energy gap vanishes at k = (7r,7r), according to 
spin-wave theory, corresponding to the expected Gold- 



stone mode. The series extrapolations in A by means of 
standard Pade approximants still give a finite result at 
that point, because they assume the series is regular at 
A = 1. The energy gap at k = (0, 0) is indeed finite, and 
behaves like yD at small \D\, mirroring that in the easy- 
axis region. It rises rapidly at large D. The qualitative 
behaviour is well reproduced by SWl theory. 



C. QPM Phase: Bulk Ground State Properties 

To investigate the large \D\ quantum paramagnetic 
phase we decompose the full Hamiltonian a,s H = Hq + 
XV, with 



-0.2 



Ho^\D\j2isn' + jJ2^^^j 



<y> 



and 



V 






Wi ^i + ^i ^j ) 



(11) 



(12) 



The unperturbed ground state has all spins in the S^ = 
state, and the effect of the perturbation is to create 
(+ -) states on neighbouring sites (exciton-antiexciton 
pairs). Note that, unlike the previous sections, we do not 
perform a spin rotation on one sublattice. The derivation 
then follows standard lines, and we have obtained series 
to order A^^ for both the ground state energy and for the 
'quadrupole moment' Q =< 3(5'f )^ — 2 >. An alternative 
approach, in which only the anisotropy term is included 
in Hq and the full exchange term as V, is also possible. 
The expansion parameter is then J/ D. We have carried 
this through but it seems not to yield any improvement 
in precision. 

In Fig. [TT]we plot the ground state energy versus D/J 
in the QPM phase. We also include some of the data 
from the PAFM expansion (Section IIII Al above) near 
the crossover region. As can be seen, the two curves 
meet smoothly, but the precision is inadequate to distin- 
guish between a second order transition or a weak first or- 
der transition (with a small discontinuity in slope). The 
quadrupole moment Q increases smoothly from a value 
—2 at large |Z3| to approximately —1.6 at D = —5, but 
shows no interesting behaviour. 



D. QPM Phase: Excitations 

The low energy excitations, termed 'excitons' and 
'antiexcitons', arise from exciting one of the S"^ = sites 
to S^ = ±1. Such a local excitation can then propagate 
through the lattice as a well defined quasiparticle with 
energy e(k). 

Using the Hamiltonian decomposition pill2p and the 
usual linked cluster methods, we have computed series 
for the excitation energy to order A-'^". Figure [T^ shows 
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FIG. 11: Ground state energy per site at negative D/J (set- 
ting J = 1). The open circles are estimates in the QPM phase, 
and the filled squares are PAFM estimates. 



a plot along symmetry lines in the Brillouin zone for two 
values of D/J, viz. D/J — —10.0, —6.0. As is apparent, 
the energy gap at (tt, tt) is closing as D increases, and 
we expect it to vanish at the quantum phase transition 
point Dc- We found that the best way of locating the 
phase transition from the QPM to the PAFM phase was 
to perform a Dlog Fade analysis of the series in A for 
the energy gap in the QPM phase at k = (tt, tt), looking 
for the zero point. Hence we estimate the critical point, 
where the energy gap goes to zero at the physical value 
A = 1, lies at {D/J)c = —5.61(5). A similar analysis of 
the magnetization in the PAFM phase gives a somewhat 
less reliable estimate, {D/J)c = —5.7(2), which is com- 
patible with the figure above. We note that the coupled 
cluster calculation p3| gives {D/J)c = -6.97, -6.38 at 
successive levels. 

The energy gap in the QPM at momentum k = (tt, tt) 
was again estimated by forming Padc approximants in 
the variable 5 — 1— {1 — X/\cY , where Ac and a are the 
location and critical index, respectively, of the energy 
gap as a function of A, estimated by the usual Dlog Fade 
methods. The index a appeared consistently as 0.70(2). 

Figure [12] shows the resulting estimates of the energy 
gap in the QPM at this momentum, as a function of 
D / J , with a fit near the critical region of the form e(k) ex 
(5.61-D/ J)'', where the fit gives v = 0.73(3). One would 
naturally conclude that the critical indices v and a are 
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FIG. 12: Single particle dispersion in the QPM phase along 
symmetry lines in the Brillouin zone, for D/J = —10.0 (cir- 
cles) and D/J = —6.0 (squares). 



FIG. 13: The energy gap e(k) at momentum k = {n, n) in 
the QPM phase, as a function of D (J = 1). The dashed line 
is a fit to the data in the region —7.0 < D/J < —5.6, in the 
neighbourhood of the critical point. 



identical. 

A similar fit of the form Ms oc (5.61 + D/J)'^ to the 
magnetization in the PAFM phase is shown in Figure 
[HI The fit over the range —5.6 < D/J < —4.5 gives 
/? = 0.25(3), and again one would conclude that the 
magnetization indices in the variables A and D/J are 
the same. 

These indices should be compared with the expected 
values for the universality class corresponding to this 
quantum phase transition, namely those of the classical 
0(2) model in three dimensions, which are v = 0.671, 
/3 = 0.346 [2^1 . The agreement is not very good, but 
this is perhaps not surprising in view of the crude and 
indirect methods used in our estimates. 

An analytic theory for the QPM phase has been 
proposed by Papanicolaou j29|], based on a generalized 
Holstein-Primakoff transformation. This gives 



e(k) = ^D{D + 8 J7k) 
which yields for the (tt, tt) gap 



e(7r,7r) = y/D{D + 8 J) 



(13) 



(14) 



i.e. {D/J)c — —8.0, u — 0.5. These do not agree well 
with the series estimates. 



IV. CONCLUSIONS 

Magnetic systems with S = \ and with easy-axis or 
easy-plane crystal field anisotropy have become of inter- 
est again, as a result of new materials and suggestions 
of novel quantum phases. Early theoretical approaches, 
based on mean field. Green's function, and spin-wave 
approximations are of uncertain and dubious validity. 
Present day (numerical) approaches such as Quantum 
Monte Carlo and series methods allow rather precise cal- 
culation of ground state properties and of the spectrum 
of elementary excitations, and hence of energy gaps. 

We have used comprehensive series methods, for the 
first time, to study the model on the square lattice. Three 
distinct phases are identified, in agreement with previous 
work. We have compared our results with those of first 
and second order spin-wave theory. In the Ising antifer- 
romagnetic (lAFM) phase, the first order theory deviates 
substantially from the series results, but the second or- 
der theory (and even a 'modified' first order theory) are 
in quantitative agreement with scries, within one or two 
percent. In the planar antiferromagnetic (PAFM) phase, 
only a first order theory is available, which gives a rea- 
sonable description at small negative D, but fails in the 
neighbourhood of the transition to the quantum param- 
agnetic (QPM) phase. 



The transition between the planar antiferromagnetic 
and quantum paramagnetic phase is located at D / J — 
—5.61(5). The transition appears to be of second order, 
with critical indices in qualitative but not quantitative 
agreement with those of the classical 0(2) model in three 
spatial dimensions. We would not claim any contradic- 
tion here as the error bars on our estimates are rather 
large. 

The series approach followed in this paper can, of 
course, be applied equally well to other lattices. Indeed 
there is considerable interest in the one-dimensional case, 
and work on this is in progress. 



J 



Acknowledgments 



This work forms part of a research project supported 
by a grant from the Australian Research Council. We are 
grateful for the computing resources provided by the Aus- 
tralian Partnership for Advanced Computing (APAC) 
National Facility and by the Australian Centre for Ad- 
vanced Computing and Communications (ACS). 



APPENDIX A: SPIN WAVE THEORY FOR THE EASY-AXIS (D > 0) CASE 

The spin- wave approximation is well known. Nevertheless, for completeness, we give a brief summary of the second 
order theory for this model, following closely the treatment of ref. |3Ql] . 
The initial Hamiltonian 



H = jY, \StS^ + \{StS- + SiS+)] -Dj2iSD' 



(Al) 



<y> 



is expressed in terms of boson operators a^, bj on the respective sublattices via a Dyson-Maleev transformation 



(A2) 



A: 5f = 5 - ala^, S+ = V2S{1 - a|a,/25)a,, Sr ^ V^Sa} 
B: Sl^b]bj-S, S+ ^V2Sb]{l-b]bj/2S), SJ = V2Sb.j 

followed by a transformation to k-space (Bloch) operators 

-^ - \/|E ^""■"'«>^' ^. - \/|E ^""^^-^ (A3) 

k k 

where the sum is over N/2 points in the reduced Brillouin zone, giving 

H = -]-NS'^{zJ + 2D) + [S{zJ + 2D)-D]^{alai, + blb]^) + zJS^-^]^{albl + aijD^) 

k k 

T7 E A (ki - k2 - ka + k4 ) [27k3 -k4 aLi ^k^ bl^ b^^ + 71,4 al^ Ok^ Oka bk^ + lk^ al^ bl^ bl^ b^^ ] 



N 
2D 



ki,k2,k3,k4 



Y^ 5(ki + ka - ks - ^4:)[al^al^ak^ak^ + bl^bl^bi^^bk^] 

ki,k2,k3,k4 

Here z is the coordination number of the lattice and 7^ is the usual 

7k; — 1/z y^ exp ik ■ 6 — -(cos kx + cos ky) 



(A4) 



(A5) 



for the SQ lattice. 

The reader's attention is drawn to the factor [S{zJ + 2D) — D] associated with the diagonal quadratic terms. If we 
consider successive orders of spin wave theory as corresponding to decreasing powers of S, then the first-order theory 
(SWl) will retain only S{zJ + 2D). However if we include the complete term, for 5 = 1, we have {zJ + D). This is 
the origin of the modified first-order theory (SWla) discussed in Section Hi Al 
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To diagonalize the quadratic part of the Hamiltonian we use a standard Bogohubov transformation 

Ok = WkAk - V]^bI 

bk = -VkAl + u^Bk (A6) 

with Uk — cosh^kji'k — sinh^k- 
This gives, after some algebra, 

H = iV^o+^^k(4A + S^Sk)+^Fk(4B^ + AkBk) 

k k 

+ ^ E '5(ki+k2-k3-k4)F4(kl,k2,k3,k4)(i?i^<i?k3Bk4+4,4.^k3AkJ 

ki ,k2:k3,k4 

+ further normal ordered quartic terms (A7) 



where 



Eo = -{2J + D) + {AJ + D)R2~4JR3-2J{R2-R3)'^-2DRl (A8) 

fly, = [4J(l-i?2 + ^3) + £'(l-4i?2)]cosh26'k-4J(l-i?2 + i?3)7kSinh26'k (A9) 

Vk = 4J(1 - i?2 + i?3)7kCosh26lk - [4J(1 - i?2 + i^s) + -0(1 - 4i?2)] sinh26'k (AlO) 

F4(kl,k2,k3,k4) = 4J[(UkiMk2"k3«k4 +WkiWk2«k3'"k4)7k4 - 2UkiWk2'"k3"k4 7k2-k4] 

-2i:)[UkiUk2Wk3"k4 +WkiWk2Wk3'"k4] (All) 



and 



k k 

X] TkWkl'k = T7 E "^"^ ^^^^ ^^k (A12) 



k k 



and we have set z = A, S = 1. 

We may choose our parameter 6'k so that Vk = 0, i.e. 

tanh 2.k - 4J(l-i^2 + i?3)7k ^ 

4J(l-i?2 + -R3)+-D(l-4i?2) ^ ' 

Dropping the quartic terms in (|A7P then yields the second order spin-wave Hamihonian 

H = NEo + Y,(i^i^U^ + BlBu) (A14) 

k 

with 

ei = [4J(1 - i?2 + i?3) + -0(1 - 4i?2)]' - [4J(1 - i?2 + -R3)7k]' (A15) 

The magnetization is 

M = S- < ajai >= • • ■ = 1 - i?2 (A16) 

These equations can then be solved numerically. Note that the expressions (|A12|) for i?2 and R^ themselves involve 
i?2 and i?3 on the right-hand side, and must be solved iteratively. A convenient starting point is the first-order 
spin- wave results. We used a double Gaussian quadrature procedure to carry out the Brillouin zone integrations. 

We can obtain the qualitative behaviour of these quantities at small D from the SWla approximation. Then 



ek - V(4J + -D)2-(4J7k)2 

~ 2V2JD as D-^0 (A17) 
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^' 8^2/^ J^ '^^^'^^^[(l + i?/4J)2-(cosA:.+cosfc,)V4] 2 

~ 0.1966- - — OS D^O (A18) 

TT \2J J 



for k — (0, 0), showing that the energy gap behaves hke \/D at small D. In the same approximation, we find 

1 f'- r^- ,^ _,,^ (i + i?/4j) 1 



using the results of [30| . Thus we see a vD singularity emerging in the magnetization near D ^ as well. 
At large D, on the other hand, we have 

Uk ~ 1, Wk ~ as D ^ 00 (A19) 

and hence in leading order the single-magnon energy is 

e^^ D + U as D -> 00 (A20) 

The 2-particle transition amplitude is 

y4(ki,k2,k3,k4) 2D as Z? -> 00 (A21) 

In terms of the centre-of-niass and relative momenta 

K = ki + k2 = kg + k4 (A22) 

q=^(ki-k2), p=i(k3-k4) (A23) 

the 2-particle bound state obeys the integral Bethe-Salpeter equation 

[E2 (K) - E, (K/2 + q) - El (K/2 - q)] V2(K, q) = ^ ^ Af (K, q, p)V'2 (K, p) (A24) 

p 

where 

M(K,q,p) = 2T/4(K,q,p) AD as D ^ 00. (A25) 

This equation is satisfied by a solution where ^2(K,q) is independent of q (corresponding to two particle excitations 
at the same point), with 

E2{K) = Ei(K/2 + q)+Ei(K/2-q)-2D 

^2{4:J + D)-2D^8J as D ^ 00 (A26) 

This is precisely the energy one would naively expect for a AS^ = ±2 excitation in this limit. 

APPENDIX B: SPIN WAVE THEORY FOR THE PAFM PHASE 

To derive spin- wave theories for the easy-plane small \D\ phase, we assume 2-sublattice Neel order in the z-direction 
and write the Hamiltonian as 

H^jY, [S^S^ + ^{S+Sr + SrS+)] IdJ2{S+ + S-f (Bl) 

where we have chosen the x-axis to be the hard direction. 

To leading order in S we again introduce boson operators in the respective sublattices 



Sf = S ~ ajoi, 5*,+ = V^ai, S, = V2Sa. 



t!,. C C+ — , /oe^t C- _ . /OCJ, V ) 



B: Sj = b]bj-S, SJ = V2Sb], SJ = V2Sbj 



12 
followed by a transformation to Bloch operators as before. Keeping only quadratic terms yields 
H = -(2J+ii?)iV + (4J-7^)^(aJ,ak + 6i,6k)+4J^7k(al,&l, + ak6k 

k k 

-2^X!("k"-k + aka-k + fo[&^k + ^k&-k (B3) 

k 

where we have set 5 = 1, z = 4. The Brillouin zone sums are over N/2 points in the reduced zone. 

To diagonalize this, more general, quadratic Hamiltonian we introduce operators {gi, ^2, 93, 94} = {flki ^ki "Lk' ^-k} 
and write the Hamiltonian as 

H^-UN+\Y,h,,4q, (B4) 



2 

k 



where h is the 4x4 matrix 



/ AJ - D 4J7k -D \ 

4J7k U^D -D 

^ \ D U~D 4 J7k ^^^^ 

V -Z? 4J7k 4J-D/ 

We now introduce a transformation to new coordinates {Qt : i = 1,4} 

q, = Y.S,,Qj (B6) 

i 

with the constraint [Qi, Qj] — Ji^ij with J^ = (1, —1, —1, 1), i.e. SJS'^ = J (where J is a diagonal matrix with entries 
•hi = Ji)- Following the argument of Tsallis [31], one easily shows that the matrix 

/ U - D -4 J7k D \ 

h-hT- 4J7k -{U-D) -D 

^ ^^ \ ^D -{U-D) 4J7k ^^^' 

V Z? -4J7k AJ-d) 

can be diagonalized by a similarity transformation, with its eigenvalues remaining invariant. It has eigenvalues 
±Ai, ±A2 where Ai, A2 are the spin- wave energies. 
The diagonalized Hamiltonian can then be written as 

H^NE^+Y, [Aik4^k + A2kSt Bk] (B8) 



with 



Direct calculation gives 



^o = -4J+-^^(Aik + A2k) (B9) 



2N 

k 



A?k - 16j2(l-72)-8i?J(l + 7k) (BIO) 

AL = 16j2(l-7^)-8i?J(l-7k) (Bll) 

In the reduced zone we have two branches, one of which is gapless at k = (0, 0) and the other at (tt, tt). However we 
note that A2(7r — fcj;,7r — ky) = \i{kx, ky) and hence in a full zone we need only consider a single branch Wk — Aik. 
Then we find that the spin wave energy vanishes at k = (tt, tt), corresponding to the expected Goldstonc mode, while 
at k = (0, 0) the gap is 4-y/J(— D), mirroring the square root behaviour found in the easy-axis case (modulo the factor 
rj referred to previously) . 

The magnetization is given by 

Af = l-l^[|^i2(k)P + |5i3(k)p] (B12) 

k 
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and can be obtained numerically from the transformation equations. 

The theory described above follows from either the Holstein-Primakoff or Dyson-Maleev approach, at lowest order. 
However an attempt to extend these to higher order fails, as the resulting spin wave energies do not possess the 
Goldstone mode required by symmetry. 
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